#### prepare packages and functions ####
require(estimatr)
require(Rcpp)
require(RcppArmadillo)
require(coda)
require(dichromat)

issue.var.name <- c("a.article9", "b.defense", "c.revisionism", "d.women", 
                    "e.gay", "f.immigrant", "g.growth", "h.tax")

color <- colorschemes$Categorical.12
plot.col <- c(color[12], color[12], color[12], 
              color[10], color[10], color[10], 
              color[2], color[2])

blue.orange <- colorschemes$BluetoOrange.10
group.col.1 <- c(blue.orange[10], "gray30", blue.orange[1])
group.col.2 <- c(paste0(blue.orange[7], "80"), "#CCCCCC80", 
                 paste0(blue.orange[4], "80"))
group.lty <- c(1, 2, 4)

issue.label <- c("Revision of Article 9", "Increasing Defense Power", 
                 "Historical Revisionism", "Women's Empowerment", 
                 "Gay Marriage", "Accepting Foreign Workers", 
                 "Economic Growth", "Progressive Tax")

# function to compute MMs
MM <- function(data, attribute, outcome, id) {
  lm.data <- data.frame(id = data[, id], A = data[, attribute], 
                        Y = data[, outcome])
  l <- nlevels(lm.data$A)
  result.matrix <- matrix(NA, l, 3)
  for (i in 1:l) {
    level.label <- levels(lm.data$A)[i]
    result.matrix[i, ] <- unlist(lm_robust(Y ~ 1, data = lm.data, 
                                           subset = A == level.label, 
                                           clusters = id)[c(1, 6, 7)])
  }
  result.matrix
}

# function for the quadratic regression
MM.predict <- function(data, attribute, outcome, predictor, values, id) {
  lm.data <- data.frame(id = data[, id], A = data[, attribute],
                        Y = data[, outcome], X = data[, predictor])
  l <- nlevels(lm.data$A)
  N <- max(lm.data$id)
  lm.result <- array(NA, c(3, 3, l))
  prediction <- array(NA, c(length(values), l, 3))
  for (i in 1:l) {
    level.label <- levels(lm.data$A)[i]
    QR.result <- lm_robust(Y ~ X + I(X ^ 2), data = lm.data,
                           subset = A == level.label,
                           clusters = id)
    lm.result[, 1, i] <- QR.result$coefficients
    lm.result[, 2, i] <- QR.result$std.error
    lm.result[, 3, i] <- QR.result$p.value
    prediction[, i, ] <- predict(QR.result, data.frame(X = values),
                                 interval = "confidence")$fit
  }
  list(lm.result = lm.result, prediction = prediction)
}

# function to conduct F-tests for linear models
age.1.F.test <- function(variable, data, choice) {
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ age + ", variable))
  } else {
    lm.formula <- as.formula(paste0("rating ~ age + ", variable))
  }
  lm.model.1 <- lm(lm.formula, data = data)
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ age * ", variable))
  } else {
    lm.formula <- as.formula(paste0("rating ~ age * ", variable))
  }
  lm.model.2 <- lm(lm.formula, data = data)
  anova.result <- anova(lm.model.1, lm.model.2)
  c(anova.result$"F"[2], anova.result$"Pr(>F)"[2])
}

# function to conduct F-tests for quadratic models
age.2.F.test <- function(variable, data, choice) {
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ age + I(age ^ 2) + ", variable))
  } else {
    lm.formula <- as.formula(paste0("rating ~ age + I(age ^ 2) + ", variable))
  }
  lm.model.1 <- lm(lm.formula, data = data)
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ (age + I(age ^ 2)) * ", variable))
  } else {
    lm.formula <- as.formula(paste0("rating ~ (age + I(age ^ 2)) * ", variable))
  }
  lm.model.2 <- lm(lm.formula, data = data)
  anova.result <- anova(lm.model.1, lm.model.2)
  c(anova.result$"F"[2], anova.result$"Pr(>F)"[2])
}

# function to estimate the finite mixture model of linear regressions
sourceCpp("finite_mixture_linear_regressions.cpp")

# function to deal with label switching
labeling <- function(mcmc.result) {
  MCMCsample <- list()
  for (i in 1:length(mcmc.result)) {
    MCMCsample[[i]] <- mcmc.result[[i]]$beta_r[1, , ]
  }
  l <- length(MCMCsample)
  g <- ncol(mcmc.result[[1]]$beta_r)
  label <- matrix(NA, l, g)
  for (i in 1:l) {
    intercept <- rowMeans(MCMCsample[[i]])
    for (j in 1:g) {
      label[i, j] <- which(c(g - rank(intercept) + 1) == j)
    }
  }
  label
}

# function to create mcmc.list objects from an Rcpp output
mcmc.list.create <- function(mcmc.result, parameter, label = NULL, group = NULL) {
  combined.list <- mcmc.list()
  if (parameter %in% c("beta_r", "beta_c")) {
    if (is.null(label)) {
      for (i in 1:length(mcmc.result)) {
        combined.list[[i]] <- mcmc(t(mcmc.result[[i]][[parameter]]))
      }
    } else {
      for (i in 1:length(mcmc.result)) {
        combined.list[[i]] <- mcmc(t(mcmc.result[[i]][[parameter]][, label[i, group], ]))
      }
    }
  } else {
    if (is.null(label)) {
      for (i in 1:length(mcmc.result)) {
        combined.list[[i]] <- mcmc(mcmc.result[[i]][[parameter]][, ])
      }
    } else {
      for (i in 1:length(mcmc.result)) {
        combined.list[[i]] <- mcmc(mcmc.result[[i]][[parameter]][, label[i, group]])
      }
    }
  }
  combined.list 
}

# function to compute the WAIC
WAIC <- function(loglik) {
  lppd <- sum(log(colMeans(exp(loglik))))
  p.WAIC2 <- sum(apply(loglik, 2, var))
  -2 * (lppd - p.WAIC2)
}

# function to summarize posterior distributions
posterior.summary <- function(x, CI = TRUE, prob = 0.95) {
  posterior.mean <- mean(x)
  if (CI == TRUE) {
    HPD.interval <- HPDinterval(as.mcmc(x), prob = prob)
    c(posterior.mean, HPD.interval)
  } else {
    posterior.mean
  }
}

#### load data ####
# task-level dataset including satisficing responses
task.raw.data <- read.csv("task_raw_data.csv")

# task-level dataset excluding satisficing responses
task.data <- read.csv("task_data.csv")

# reorder the levels of attribute variables
position.label <- c("Agree", "Neither", "Disagree")
task.raw.data$a.article9 <- factor(task.raw.data$a.article9, levels = position.label)
task.raw.data$b.defense <- factor(task.raw.data$b.defense, levels = position.label)
task.raw.data$c.revisionism <- factor(task.raw.data$c.revisionism, levels = position.label)
task.raw.data$d.women <- factor(task.raw.data$d.women, levels = position.label)
task.raw.data$e.gay <- factor(task.raw.data$e.gay, levels = position.label)
task.raw.data$f.immigrant <- factor(task.raw.data$f.immigrant, levels = position.label)
task.raw.data$g.growth <- factor(task.raw.data$g.growth, levels = position.label)
task.raw.data$h.tax <- factor(task.raw.data$h.tax, levels = position.label)
task.data$a.article9 <- factor(task.data$a.article9, levels = position.label)
task.data$b.defense <- factor(task.data$b.defense, levels = position.label)
task.data$c.revisionism <- factor(task.data$c.revisionism, levels = position.label)
task.data$d.women <- factor(task.data$d.women, levels = position.label)
task.data$e.gay <- factor(task.data$e.gay, levels = position.label)
task.data$f.immigrant <- factor(task.data$f.immigrant, levels = position.label)
task.data$g.growth <- factor(task.data$g.growth, levels = position.label)
task.data$h.tax <- factor(task.data$h.tax, levels = position.label)

original.data <- subset(read.csv("ideological_label_conjoint_data.csv", 
                                 fileEncoding = "utf8"), 
                        Q5.2 < 53)
N <- nrow(original.data)  # number of respondents including satisficers

#### general understanding ####
## compute MMs
# including satisficing responses
choice.result.entire <- rating.result.entire <- matrix(NA, 24, 3)
for (i in 1:8) {
  choice.result.entire[(3 * (i - 1) + 1):(3 * i), ] <- MM(task.raw.data, 
                                                          issue.var.name[i], 
                                                          "choice", 
                                                          "respondent.id")
  rating.result.entire[(3 * (i - 1) + 1):(3 * i), ] <- MM(task.raw.data, 
                                                          issue.var.name[i], 
                                                          "rating", 
                                                          "respondent.id")
}

# excluding satisficing responses
choice.result <- rating.result <- matrix(NA, 24, 3)
for (i in 1:8) {
  choice.result[(3 * (i - 1) + 1):(3 * i), ] <- MM(task.data, 
                                                   issue.var.name[i], 
                                                   "choice", 
                                                   "respondent.id")
  rating.result[(3 * (i - 1) + 1):(3 * i), ] <- MM(task.data, 
                                                   issue.var.name[i], 
                                                   "rating", 
                                                   "respondent.id")
}

# Figure A.5
cairo_pdf("Figure_A5.pdf", width = 6.2, height = 4, pointsize = 8)
par(mar = c(3, 0, 1.5, 0), lwd = 0.5, xpd = TRUE)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.23, 0.85), ylim = c(0, 33), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(seq(0.4, 0.6, 0.05), 0, seq(0.4, 0.6, 0.05), 33, lwd = 0.4, col = "gray")
for (i in 1:8) {
  segments(0.395, 35 - 4 * i, 0.605, 35 - 4 * i, lty = 3, col = "gray")
  segments(0.395, 34 - 4 * i, 0.605, 34 - 4 * i, lty = 3, col = "gray")
  segments(0.395, 33 - 4 * i, 0.605, 33 - 4 * i, lty = 3, col = "gray")
  segments(choice.result.entire[(3 * (i - 1) + 1):(3 * i), 2], 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) + 0.2, 
           choice.result.entire[(3 * (i - 1) + 1):(3 * i), 3], 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) + 0.2, col = plot.col[i])
  points(choice.result.entire[(3 * (i - 1) + 1):(3 * i), 1], 
         c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) + 0.2, 
         pch = 19, col = plot.col[i])
  segments(choice.result[(3 * (i - 1) + 1):(3 * i), 2], 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) - 0.2, 
           choice.result[(3 * (i - 1) + 1):(3 * i), 3], 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) - 0.2, col = plot.col[i])
  points(choice.result[(3 * (i - 1) + 1):(3 * i), 1], 
         c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) - 0.2, 
         pch = 21, col = plot.col[i], bg = "white")
}
axis(1, at = seq(0.4, 0.6, 0.05))
segments(seq(0.65, 0.85, 0.05), 0, seq(0.65, 0.85, 0.05), 33, lwd = 0.4, col = "gray")
for (i in 1:8) {
  segments(0.645, 35 - 4 * i, 0.855, 35 - 4 * i, lty = 3, col = "gray")
  segments(0.645, 34 - 4 * i, 0.855, 34 - 4 * i, lty = 3, col = "gray")
  segments(0.645, 33 - 4 * i, 0.855, 33 - 4 * i, lty = 3, col = "gray")
  segments(0.75 + (rating.result.entire[(3 * (i - 1) + 1):(3 * i), 2] - 5.5) / 4, 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) + 0.2, 
           0.75 + (rating.result.entire[(3 * (i - 1) + 1):(3 * i), 3] - 5.5) / 4, 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) + 0.2, col = plot.col[i])
  points(0.75 + (rating.result.entire[(3 * (i - 1) + 1):(3 * i), 1] - 5.5) / 4, 
         c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) + 0.2, 
         pch = 19, col = plot.col[i])
  segments(0.75 + (rating.result[(3 * (i - 1) + 1):(3 * i), 2] - 5.5) / 4, 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) - 0.2, 
           0.75 + (rating.result[(3 * (i - 1) + 1):(3 * i), 3] - 5.5) / 4, 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) - 0.2, col = plot.col[i])
  points(0.75 + (rating.result[(3 * (i - 1) + 1):(3 * i), 1] - 5.5) / 4, 
         c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i) - 0.2, 
         pch = 21, col = plot.col[i], bg = "white")
}
axis(1, at = seq(0.65, 0.85, 0.05), labels = seq(5.1, 5.9, 0.2))
for (i in 1:8) {
  text(0.23, 36 - 4 * i, labels = issue.label[i], pos = 4)
  text(0.33, 35 - 4 * i, labels = "Agree", cex = 0.9, pos = 4)
  text(0.33, 34 - 4 * i, labels = "Neither", cex = 0.9, pos = 4)
  text(0.33, 33 - 4 * i, labels = "Disagree", cex = 0.9, pos = 4)
}
legend(0.23, 0, c("w/ satisficing responses", "w/o satisficing responses"), 
       pch = c(19, 21), bg = c(NA, "white"), bty = "n")
mtext(c("Choice", "Rating"), at = c(0.5, 0.75), cex = 1.2, font = 2)
dev.off()

#### understanding conditioned by age ####
## predict MMs by age (divided by 10)
choice.predict.entire <- rating.predict.entire <- list()
for (i in 1:8) {
  choice.predict.entire[[i]] <- MM.predict(task.raw.data, issue.var.name[i], 
                                           "choice","age.decimal", 
                                           c(18:69) / 10, "respondent.id")
  rating.predict.entire[[i]] <- MM.predict(task.raw.data, issue.var.name[i], 
                                           "rating", "age.decimal", 
                                           c(18:69) / 10, "respondent.id")
}

# Figure A.6(a)
cairo_pdf("Figure_A6a.pdf", width = 6.2, height = 3.3, pointsize = 8)
layout(matrix(1:8, 2, 4, byrow = TRUE))
par(mar = c(3, 3, 3, 1), lwd = 0.5)
for (i in 1:8) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(18, 69), ylim = c(0.3, 0.7), 
       main = issue.label[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  for (j in 3:1) {
    polygon(c(18:69, 69:18), 
            c(choice.predict.entire[[i]]$prediction[, j, 2], 
              rev(choice.predict.entire[[i]]$prediction[, j, 3])), 
            border = NA, col = group.col.2[j])
  }
  for (j in 3:1) {
    lines(18:69, choice.predict.entire[[i]]$prediction[, j, 1], 
          lty = j, col = group.col.1[j])
  }
  if (i == 1) {
    text(69, 0.61, "Agree", pos = 2)
    text(69, 0.51, "Neither", pos = 2)
    text(69, 0.38, "Disagree", pos = 2)
  }
  axis(1, at = c(seq(18, 58, 10), 69), lwd = 0.5)
  axis(2, at = seq(0.3, 0.7, 0.1), lwd = 0.5)
  mtext("Age", side = 1, line = 2, cex = 0.7)
  mtext("Marginal Mean", side = 2, line = 2, cex = 0.7)
}
dev.off()

# Figure A.6(b)
cairo_pdf("Figure_A6b.pdf", width = 6.2, height = 3.3, pointsize = 8)
layout(matrix(1:8, 2, 4, byrow = TRUE))
par(mar = c(3, 3, 3, 1), lwd = 0.5)
for (i in 1:8) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(18, 69), ylim = c(4.7, 6.3), 
       main = issue.label[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  for (j in 3:1) {
    polygon(c(18:69, 69:18), 
            c(rating.predict.entire[[i]]$prediction[, j, 2], 
              rev(rating.predict.entire[[i]]$prediction[, j, 3])), 
            border = NA, col = group.col.2[j])
  }
  for (j in 3:1) {
    lines(18:69, rating.predict.entire[[i]]$prediction[, j, 1], 
          lty = j, col = group.col.1[j])
  }
  if (i == 1) {
    text(69, 6.27, "Agree", pos = 2)
    text(69, 5.72, "Neither", pos = 2)
    text(69, 5.15, "Disagree", pos = 2)
  }
  axis(1, at = c(seq(18, 58, 10), 69), lwd = 0.5)
  axis(2, at = seq(4.7, 6.3, 0.4), lwd = 0.5)
  mtext("Age", side = 1, line = 2, cex = 0.7)
  mtext("Marginal Mean", side = 2, line = 2, cex = 0.7)
}
dev.off()

#### latent groups ####
## prepare variables
# individual average of rating responses
average.rating.entire <- rep(NA, N)
for (i in unique(task.raw.data$respondent.id)) {
  average.rating.entire[i] <- 
    mean(task.raw.data$rating[task.raw.data$respondent.id == i])
}
task.raw.data$average.rating <- average.rating.entire[task.raw.data$respondent.id]

# dependent variables
y.c.entire <- task.raw.data$choice
y.r.entire <- task.raw.data$rating - task.raw.data$average.rating

# independent variables
lm.result.entire <- lm(rating ~ a.article9 + b.defense + c.revisionism + d.women + 
                         e.gay + f.immigrant + g.growth + h.tax, 
                       data = task.raw.data, x = TRUE)
x.entire <- lm.result.entire$x[, -1]

# rearrange respondent ids
id.entire <- as.numeric(as.factor(task.raw.data$respondent.id))

# number of respondents
l.entire <- length(unique(id.entire))

## estimation
# CAUTION: It took over 15 hours to estimate the models in the authors'
# environment (using a high performance PC as of 2021).
# If you want to skip the estimation process, request Rdata files
# containing the MCMC results from the corresponding author.
for (i in 1:6) {
  g <- i
  FMR.result.entire <- list()
  for (j in 1:3) {
    set.seed(1000 * j)
    FMR.result.entire[[j]] <- FMR(y_r = y.r.entire,  # rating-based outcomes
                                  y_c = y.c.entire,  # choice-based outcomes
                                  x = x.entire,  # independent variables
                                  id = id.entire,  # respondent ids
                                  l = l.entire,  # number of respondents
                                  g = g,  # number of latent groups
                                  b0_r = rep(0, ncol(x.entire) + 1),  # prior means for rating coefficients
                                  B0_r = diag(100 ^ 2, ncol(x.entire) + 1),  # prior variance matrix for rating coefficients
                                  b0_c = rep(0, ncol(x.entire) + 1),  # prior means for choice coefficients
                                  B0_c = diag(100 ^ 2, ncol(x.entire) + 1),  # prior variance matrix for choice coefficients
                                  v0_r = 0.01,  # prior degrees of freedom for the error variance of rating
                                  s0_r = 1,  # prior scale for the error variance of rating
                                  v0_c = 0.01,  # prior degrees of freedom for the error variance of choice
                                  s0_c = 1,  # prior scale for the error variance of choice
                                  w0 = rep(1, g),  # prior concentration parameter
                                  burnin = 5000,  # number of burn-in iterations
                                  iter = 10000,  # number of sampling iterations
                                  thin = 10  # thinning interval
    )
  }
  print(paste0("Estimation of the ", i, "-group model has been finished at ", date(), "."))
  save(FMR.result.entire, file = paste0("FMR_result_", i, "_entire.Rdata"))
  rm(FMR.result.entire)
  gc(); gc()
}

# check convergence and compute the WAICs
for (i in 1:6) {
  load(paste0("FMR_result_", i, "_entire.Rdata"))
  # relabel to deal with between-chain label switching
  if (i > 1) {
    label <- labeling(FMR.result.entire)
  } else {
    label <- matrix(1, 3, 1)
  }
  print(paste0("Convergence check for the ", i, "-group model:"))
  for (j in 1:i) {
    print(sum(gelman.diag(mcmc.list.create(FMR.result.entire, "beta_r", label, j), 
                          multivariate = FALSE)$psrf[, 1] > 1.1))
    print(sum(gelman.diag(mcmc.list.create(FMR.result.entire, "beta_c", label, j), 
                          multivariate = FALSE)$psrf[, 1] > 1.1))
    print(sum(gelman.diag(mcmc.list.create(FMR.result.entire, "sigma_r", label, j), 
                          multivariate = FALSE)$psrf[, 1] > 1.1))
    print(sum(gelman.diag(mcmc.list.create(FMR.result.entire, "sigma_c", label, j), 
                          multivariate = FALSE)$psrf[, 1] > 1.1))
  }
  print(paste0("The ", i, "-group model's WAIC:"))
  print(round(WAIC(rbind(FMR.result.entire[[1]]$loglik, 
                         FMR.result.entire[[2]]$loglik, 
                         FMR.result.entire[[3]]$loglik))))
  rm(FMR.result.entire)
  gc(); gc()
}

## post-processing
# adopt the four-group model
load("FMR_result_4_entire.Rdata")

# relabel to deal with between-chain label switching
label.entire <- labeling(FMR.result.entire)

# reorder latent groups to make interpretations easier
beta_r.sample.entire <- array(NA, c(17, 4, 30000))
for (i in 1:4) {
  beta_r.sample.entire[, i, ] <- cbind(FMR.result.entire[[1]]$beta_r[, label.entire[1, i], ], 
                                       FMR.result.entire[[2]]$beta_r[, label.entire[2, i], ], 
                                       FMR.result.entire[[3]]$beta_r[, label.entire[3, i], ])
}
round(apply(beta_r.sample.entire, c(1, 2), mean), 3)
group.order.entire <- c(1, 4, 2, 3)  # reorder groups based on the coefficients of Article 9 in absolute value

# group assignment indices
S.sample.entire <- array(NA, c(l.entire, 4, 30000))
for (i in 1:4) {
  S.sample.entire[, i, ] <- cbind(FMR.result.entire[[1]]$S[, label.entire[1, i], ], 
                                  FMR.result.entire[[2]]$S[, label.entire[2, i], ], 
                                  FMR.result.entire[[3]]$S[, label.entire[3, i], ])
}
S.sample.entire <- S.sample.entire[, group.order.entire, ]

# obtain group-wise estimates of MMs by Monte Carlo marginalization
MM.posterior.sim.choice.entire <- MM.posterior.sim.rating.entire <- 
  array(NA, c(24, 4, 30000))
for (i in 1:30000) {
  for (j in 1:4) {
    # detect observations assigned to group j in iteration i
    temp.data <- task.raw.data[id.entire %in% which(S.sample.entire[, j, i] == 1), ]
    n.obs <- nrow(temp.data)
    # posterior draw of the MM for the three levels of attribute k
    for (k in 1:8) {
      temp.mean <- tapply(temp.data$choice, temp.data[, issue.var.name[k]], mean)
      temp.sd <- tapply(temp.data$choice, temp.data[, issue.var.name[k]], sd) / sqrt(n.obs)
      MM.posterior.sim.choice.entire[(3 * (k - 1) + 1):(3 * (k - 1) + 3), j, i] <- 
        rnorm(3, temp.mean, temp.sd)
      temp.mean <- tapply(temp.data$rating, temp.data[, issue.var.name[k]], mean)
      temp.sd <- tapply(temp.data$rating, temp.data[, issue.var.name[k]], sd) / sqrt(n.obs)
      MM.posterior.sim.rating.entire[(3 * (k - 1) + 1):(3 * (k - 1) + 3), j, i] <- 
        rnorm(3, temp.mean, temp.sd)
    }
  }
}

# Figure A.7(a)
MM.posterior.summary.choice.entire <- 
  apply(MM.posterior.sim.choice.entire, 1:2, posterior.summary)

cairo_pdf("Figure_A7a.pdf", width = 6.2, height = 4.4, pointsize = 8)
par(mar = c(2, 0, 1, 0))
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.97, 3.05), ylim = c(0, 33), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
for (i in 1:4) {
  segments(seq(0.2 + 0.75 * (i - 1), 0.8 + 0.75 * (i - 1), 0.15), 0, 
           seq(0.2 + 0.75 * (i - 1), 0.8 + 0.75 * (i - 1), 0.15), 33, lwd = 0.4, col = "gray")
  for (j in 1:8) {
    segments(0.17 + 0.75 * (i - 1), 35 - 4 * j, 
             0.83 + 0.75 * (i - 1), 35 - 4 * j, lty = 3, col = "gray")
    segments(0.17 + 0.75 * (i - 1), 34 - 4 * j, 
             0.83 + 0.75 * (i - 1), 34 - 4 * j, lty = 3, col = "gray")
    segments(0.17 + 0.75 * (i - 1), 33 - 4 * j, 
             0.83 + 0.75 * (i - 1), 33 - 4 * j, lty = 3, col = "gray")
    points(MM.posterior.summary.choice.entire[1, (3 * (j - 1) + 1):(3 * j), i] + 0.75 * (i - 1), 
           c(35 - 4 * j, 34 - 4 * j, 33 - 4 * j), 
           pch = 19, col = plot.col[j])
    segments(MM.posterior.summary.choice.entire[2, (3 * (j - 1) + 1):(3 * j), i] + 0.75 * (i - 1), 
             c(35 - 4 * j, 34 - 4 * j, 33 - 4 * j), 
             MM.posterior.summary.choice.entire[3, (3 * (j - 1) + 1):(3 * j), i] + 0.75 * (i - 1), 
             c(35 - 4 * j, 34 - 4 * j, 33 - 4 * j), col = plot.col[j])
  }
  axis(1, at = seq(0.2 + 0.75 * (i - 1), 0.8 + 0.75 * (i - 1), 0.15), 
       labels = c("0.2", "", "0.5", "", "0.8"))
}
for (i in 1:8) {
  text(-0.97, 36 - 4 * i, labels = issue.label[i], pos = 4)
  text(-0.31, 35 - 4 * i, labels = "Agree", cex = 0.9, pos = 4)
  text(-0.31, 34 - 4 * i, labels = "Neither", cex = 0.9, pos = 4)
  text(-0.31, 33 - 4 * i, labels = "Disagree", cex = 0.9, pos = 4)
}
mtext(c("Group 1", "Group 2", "Group 3", "Group 4"), 
      at = seq(0.5, 2.75, 0.75), line = -0.5, cex = 1.2)
dev.off()

# Figure A.7(b)
MM.posterior.summary.rating.entire <- apply(MM.posterior.sim.rating.entire, 1:2, posterior.summary)

cairo_pdf("Figure_A7b.pdf", width = 6.2, height = 4.4, pointsize = 8)
par(mar = c(2, 0, 1, 0))
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.6, 14), ylim = c(0, 33), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
for (i in 1:4) {
  segments(seq(4.5 + 2.5 * (i - 1), 6.5 + 2.5 * (i - 1), 0.5), 0, 
           seq(4.5 + 2.5 * (i - 1), 6.5 + 2.5 * (i - 1), 0.5), 33, lwd = 0.4, col = "gray")
  for (j in 1:8) {
    segments(4.4 + 2.5 * (i - 1), 35 - 4 * j, 
             6.6 + 2.5 * (i - 1), 35 - 4 * j, lty = 3, col = "gray")
    segments(4.4 + 2.5 * (i - 1), 34 - 4 * j, 
             6.6 + 2.5 * (i - 1), 34 - 4 * j, lty = 3, col = "gray")
    segments(4.4 + 2.5 * (i - 1), 33 - 4 * j, 
             6.6 + 2.5 * (i - 1), 33 - 4 * j, lty = 3, col = "gray")
    points(MM.posterior.summary.rating.entire[1, (3 * (j - 1) + 1):(3 * j), i] + 2.5 * (i - 1), 
           c(35 - 4 * j, 34 - 4 * j, 33 - 4 * j), 
           pch = 19, col = plot.col[j])
    segments(MM.posterior.summary.rating.entire[2, (3 * (j - 1) + 1):(3 * j), i] + 2.5 * (i - 1), 
             c(35 - 4 * j, 34 - 4 * j, 33 - 4 * j), 
             MM.posterior.summary.rating.entire[3, (3 * (j - 1) + 1):(3 * j), i] + 2.5 * (i - 1), 
             c(35 - 4 * j, 34 - 4 * j, 33 - 4 * j), col = plot.col[j])
  }
  axis(1, at = seq(4.5 + 2.5 * (i - 1), 6.5 + 2.5 * (i - 1), 0.5), 
       labels = c("4.5", "", "5.5", "", "6.5"))
}
for (i in 1:8) {
  text(0.6, 36 - 4 * i, labels = issue.label[i], pos = 4)
  text(2.8, 35 - 4 * i, labels = "Agree", cex = 0.9, pos = 4)
  text(2.8, 34 - 4 * i, labels = "Neither", cex = 0.9, pos = 4)
  text(2.8, 33 - 4 * i, labels = "Disagree", cex = 0.9, pos = 4)
}
mtext(c("Group 1", "Group 2", "Group 3", "Group 4"), 
      at = seq(5.5, 13, 2.5), line = -0.5, cex = 1.2)
dev.off()

# estimation results of parameters other than coefficients (Table A.8)
# error variance for choice-based outcomes
sigma_c.sample.entire <- matrix(NA, 30000, 4)
for (i in 1:4) {
  sigma_c.sample.entire[, i] <- rbind(FMR.result.entire[[1]]$sigma_c[, label.entire[1, i]], 
                                      FMR.result.entire[[2]]$sigma_c[, label.entire[2, i]], 
                                      FMR.result.entire[[3]]$sigma_c[, label.entire[3, i]])
}
sigma_c.sample.entire <- sigma_c.sample.entire[, group.order.entire]
sigma_c.summary.entire <- apply(sigma_c.sample.entire, 2, posterior.summary, CI = TRUE)
round(sigma_c.summary.entire, 2)

# error variance for rating-based outcomes
sigma_r.sample.entire <- matrix(NA, 30000, 4)
for (i in 1:4) {
  sigma_r.sample.entire[, i] <- rbind(FMR.result.entire[[1]]$sigma_r[, label.entire[1, i]], 
                                      FMR.result.entire[[2]]$sigma_r[, label.entire[2, i]], 
                                      FMR.result.entire[[3]]$sigma_r[, label.entire[3, i]])
}
sigma_r.sample.entire <- sigma_r.sample.entire[, group.order.entire]
sigma_r.summary <- apply(sigma_r.sample.entire, 2, posterior.summary, CI = TRUE)
round(sigma_r.summary, 2)

# weight parameters
eta.sample.entire <- matrix(NA, 30000, 4)
for (i in 1:4) {
  eta.sample.entire[, i] <- rbind(FMR.result.entire[[1]]$eta[, label.entire[1, i]], 
                                  FMR.result.entire[[2]]$eta[, label.entire[2, i]], 
                                  FMR.result.entire[[3]]$eta[, label.entire[3, i]])
}
eta.sample.entire <- eta.sample.entire[, group.order.entire]
eta.summary <- apply(eta.sample.entire, 2, posterior.summary, CI = TRUE)
round(eta.summary, 2)

# proportion of latent groups
round(apply(colMeans(S.sample.entire), 1, posterior.summary), 2)